home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / SoundAndMusic / cmix / lib / hplset.f < prev    next >
Encoding:
Text File  |  1988-01-31  |  5.8 KB  |  92 lines
  1.       subroutine hplset(xxlp,dur,dynam,plamp,seed,sr,q)                  
  2. c...............need to replace rand unit              
  3. c.................compile with order -lF77 -lm    
  4.       dimension q(1)                                                    
  5.       complex j                                                         
  6. c------------------decay stretch & shortening factors                   
  7. c    write(6,999)xlp,dur,dynam,plamp,seed,sr
  8. c999   format(6f10.4)
  9.       xlp=xxlp
  10.       freq=1./xlp                                                       
  11.       pi=3.141592654
  12.       w=2.*pi*freq                                                      
  13. c.....tau => decay time required                                        
  14.       tau=dur/6.91                                                      
  15.       n1=aint(sr/freq-.5)                                               
  16. 21    ga1=exp(-(n1+.5)/(tau*sr))                                        
  17.       x=1-ga1**2.                                                       
  18.       b=2.*(1.-cos(w/sr))                                               
  19.       s2=b**2.-4.*b*x                                                   
  20.       if(s2.lt.0.) goto 11                                              
  21. c------------------------                                               
  22.       s1=sqrt(s2)/(2*b)                                                 
  23.       s1=.5+s1                                                          
  24. c.....rho => loss factor                                                
  25.       rho=1.0                                                           
  26.       if(s1.lt.1) goto 10                                               
  27. c------------------------                                               
  28. 11    s1=.5                                                             
  29.       rho=exp(-1./(tau*freq))/abs(cos(w/(2*sr)))                        
  30. c------------------------                                               
  31. 10    c=1.-s1                                                           
  32.       c1=s1                                                             
  33. c----------------- tuning filter coefficient                            
  34.       xlp=sr/freq                                                       
  35. c.....pa => phase delay due to the 2 point average filter               
  36.       pa=(-1)*sr/w*atan(-c*sin(w/sr)/((1.-c)+c*cos(w/sr)))              
  37.       pa=1.-pa                                                          
  38. c.....n => buffer length                                                
  39.       n=aint(xlp-pa)                                                    
  40. c.....ga=sqrt((c1**2)+(c**2)+2.*c*c1*cos(w/sr))                         
  41. c.....tau=-(n+pa)/(sr*alog(ga))                                         
  42.       if(n.eq.n1) goto 20                                               
  43.       n1=n                                                              
  44.       goto 21                                                           
  45. c.....pc => delay required from the allpass filter                      
  46. 20    pc=(sr/freq)-n-pa                                                 
  47. c.....cc => coefficient of the allpass filter                           
  48.       cc=sin((w/sr-w/sr*pc)/2.)/sin((w/sr+w/sr*pc)/2.)                  
  49. c........                                                               
  50. c----------------- dynamic filter coefficient                           
  51.       fl=20.                                                            
  52.       fu=sr/2.                                                          
  53.       fm=sqrt(fl*fu)                                                    
  54.       rl=exp(-pi*dynam/sr)                                              
  55.       j=cmplx(0.,-1.)                                                   
  56.       gl=(1.-rl)/cabs(1.-(rl*cexp(j*2*pi*fm/sr)))                       
  57.       gl2=gl**2                                                         
  58.       w2=w/2                                                            
  59.       r1=(1.-gl2*cos(w/sr))/(1-gl2)                                     
  60.       r2=2*gl*sin(w2/sr)*sqrt(1.-gl2*(cos(w2/sr)**2))/(1.-gl2)          
  61.       r=r1+r2                                                           
  62.       if(r.ge.1) r=r1-r2                                                
  63. c-------------------                                                    
  64.       q(2)=n+20                                                         
  65.       len=q(2)                                                          
  66.       do 1 m=20,len                                                     
  67.         seed=frand(seed)
  68.     x=seed*2-1.
  69.         q(m)=plamp                                                      
  70.         if(x.lt.0.) q(m)=-plamp                                         
  71. 1     continue                                                          
  72. c--------------------                                                   
  73.     q(2)=q(2)-1
  74. c.......reset for C code which starts from 0
  75.       q(1)=q(2)                                                         
  76.       q(3)=c                                                            
  77.       q(4)=1.-c                                                         
  78.       q(9)=r                                                            
  79.       q(10)=rho                                                         
  80.       q(11)=cc                                                          
  81.       q(5)=0.                                                           
  82.       q(6)=0.                                                           
  83.       q(7)=0.                                                           
  84.       q(8)=0.                                                           
  85.       return                                                            
  86.       end                                                               
  87.       function frand(x)
  88.       n=x*1048576.
  89.       frand=float(mod(1061*n+221589,1048576))/1048576.
  90.       return
  91.       end
  92.